stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -449.879
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -33.767 0.000103114 0.000404158 0.825 0.825 28
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -33.77
## b0 6.79
## b1 1.89
## s 3.28
## m0[1] 15.48
## m0[2] 27.72
## m0[3] 9.45
## m0[4] 2.88
## m0[5] 29.60
## m0[6] 37.52
## m0[7] 19.87
## m0[8] 34.79
## m0[9] 11.39
## m0[10] 31.28
## m0[11] 29.52
## m0[12] 9.94
## m0[13] 26.71
## m0[14] 46.07
## m0[15] 8.25
## m0[16] 21.68
## m0[17] 25.69
## m0[18] 18.13
## m0[19] 41.38
## m0[20] 34.53
## y0[1] 16.95
## y0[2] 29.28
## y0[3] 15.41
## y0[4] 1.50
## y0[5] 29.19
## y0[6] 38.80
## y0[7] 14.03
## y0[8] 32.00
## y0[9] 19.81
## y0[10] 32.30
## y0[11] 30.46
## y0[12] 17.10
## y0[13] 29.38
## y0[14] 46.54
## y0[15] 4.22
## y0[16] 24.14
## y0[17] 25.40
## y0[18] 16.10
## y0[19] 45.04
## y0[20] 36.64
## m1[1] 2.88
## m1[2] 7.68
## m1[3] 12.48
## m1[4] 17.27
## m1[5] 22.07
## m1[6] 26.87
## m1[7] 31.67
## m1[8] 36.47
## m1[9] 41.27
## m1[10] 46.07
## y1[1] 3.97
## y1[2] 5.35
## y1[3] 14.07
## y1[4] 15.00
## y1[5] 25.49
## y1[6] 28.47
## y1[7] 33.38
## y1[8] 33.03
## y1[9] 39.83
## y1[10] 44.76
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -34.23 -33.84 1.36 1.07 -36.89 -32.78 1.01 469 618
## b0 6.92 7.01 1.57 1.40 4.25 9.60 1.01 211 215
## b1 1.88 1.88 0.14 0.13 1.65 2.11 1.01 294 361
## s 3.70 3.59 0.70 0.64 2.75 5.00 1.00 1158 1049
## m0[1] 15.57 15.59 1.09 0.99 13.79 17.41 1.02 229 303
## m0[2] 27.75 27.74 0.88 0.85 26.31 29.25 1.00 1297 1060
## m0[3] 9.56 9.62 1.41 1.27 7.21 11.96 1.02 211 246
## m0[4] 3.03 3.12 1.82 1.60 0.00 6.06 1.02 211 226
## m0[5] 29.62 29.60 0.93 0.89 28.15 31.20 1.00 2021 1349
## m0[6] 37.51 37.50 1.27 1.25 35.45 39.58 1.00 1789 1238
## m0[7] 19.94 19.93 0.93 0.85 18.47 21.47 1.01 303 443
## m0[8] 34.79 34.76 1.13 1.11 32.96 36.66 1.00 2220 1291
## m0[9] 11.50 11.57 1.30 1.16 9.34 13.72 1.02 213 265
## m0[10] 31.30 31.28 0.98 0.94 29.76 32.96 1.00 2337 1384
## m0[11] 29.54 29.52 0.92 0.89 28.07 31.11 1.00 2030 1349
## m0[12] 10.06 10.12 1.38 1.23 7.77 12.42 1.02 212 257
## m0[13] 26.75 26.73 0.87 0.84 25.34 28.24 1.00 1023 992
## m0[14] 46.01 45.99 1.78 1.72 43.08 48.92 1.00 818 981
## m0[15] 8.37 8.45 1.48 1.32 5.89 10.89 1.01 211 218
## m0[16] 21.75 21.73 0.88 0.84 20.37 23.22 1.01 381 517
## m0[17] 25.74 25.71 0.86 0.82 24.34 27.18 1.00 811 960
## m0[18] 18.21 18.22 0.98 0.92 16.65 19.86 1.01 260 350
## m0[19] 41.35 41.33 1.49 1.46 38.92 43.75 1.00 1131 1077
## m0[20] 34.53 34.50 1.12 1.10 32.73 36.40 1.00 2246 1291
## y0[1] 15.48 15.38 3.95 3.81 9.23 22.04 1.00 1695 1585
## y0[2] 27.91 27.80 3.89 3.58 21.65 34.44 1.00 1983 2009
## y0[3] 9.62 9.62 4.12 3.92 3.05 16.36 1.00 1113 1322
## y0[4] 3.18 3.20 4.25 3.88 -3.73 10.32 1.00 751 1086
## y0[5] 29.61 29.60 3.85 3.76 23.36 35.95 1.00 1968 1927
## y0[6] 37.48 37.59 3.95 3.69 30.97 44.04 1.00 2148 2057
## y0[7] 19.92 19.96 3.87 3.69 13.56 25.93 1.00 1269 1702
## y0[8] 34.60 34.64 3.89 3.69 27.88 40.94 1.00 1857 1933
## y0[9] 11.44 11.45 4.03 3.79 4.80 18.13 1.00 1183 1681
## y0[10] 31.27 31.30 3.97 3.89 24.83 37.70 1.00 1952 1943
## y0[11] 29.55 29.58 3.81 3.45 23.22 35.77 1.00 2152 1933
## y0[12] 10.03 10.00 4.02 3.87 3.44 16.48 1.00 1284 1763
## y0[13] 26.64 26.67 3.88 3.80 20.16 32.95 1.00 1766 1851
## y0[14] 45.98 46.07 4.16 4.10 39.10 52.82 1.00 1586 1655
## y0[15] 8.36 8.39 4.02 3.78 1.75 14.81 1.00 1091 887
## y0[16] 21.69 21.69 3.88 3.69 15.35 28.02 1.00 1749 1852
## y0[17] 25.93 25.89 3.88 3.75 19.61 31.98 1.00 2043 1837
## y0[18] 18.18 18.20 3.82 3.68 11.96 24.59 1.00 1763 1834
## y0[19] 41.34 41.22 4.07 3.87 34.76 47.93 1.00 1970 1892
## y0[20] 34.75 34.73 3.95 3.77 28.46 41.09 1.00 1973 1830
## m1[1] 3.03 3.12 1.82 1.60 0.00 6.06 1.02 211 226
## m1[2] 7.80 7.89 1.52 1.35 5.25 10.39 1.01 211 218
## m1[3] 12.58 12.63 1.24 1.11 10.52 14.70 1.02 216 282
## m1[4] 17.36 17.36 1.01 0.93 15.73 19.07 1.01 247 322
## m1[5] 22.13 22.12 0.88 0.83 20.76 23.59 1.00 405 540
## m1[6] 26.91 26.89 0.87 0.84 25.50 28.39 1.00 1063 1023
## m1[7] 31.69 31.67 0.99 0.96 30.11 33.35 1.00 2355 1448
## m1[8] 36.46 36.45 1.21 1.20 34.48 38.44 1.00 2006 1255
## m1[9] 41.24 41.22 1.48 1.45 38.82 43.63 1.00 1141 1077
## m1[10] 46.01 45.99 1.78 1.72 43.08 48.92 1.00 818 981
## y1[1] 3.00 2.97 4.15 4.01 -3.72 9.75 1.00 908 1302
## y1[2] 7.78 7.81 4.09 4.01 1.06 14.33 1.00 1116 1396
## y1[3] 12.64 12.49 3.93 3.84 6.33 19.04 1.00 1433 1501
## y1[4] 17.38 17.29 3.73 3.66 11.40 23.77 1.00 1746 1968
## y1[5] 22.13 22.08 3.86 3.79 15.73 28.55 1.00 1928 1786
## y1[6] 26.86 26.99 3.81 3.60 20.47 33.07 1.00 2065 1978
## y1[7] 31.55 31.56 3.78 3.69 25.25 37.72 1.00 1964 2014
## y1[8] 36.39 36.46 4.03 3.76 29.68 42.75 1.00 1988 1891
## y1[9] 41.43 41.48 4.01 3.74 34.60 47.90 1.00 1851 1868
## y1[10] 45.96 45.94 4.15 4.01 38.96 52.71 1.00 1899 1966
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 finished in 0.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -34.52 -36.49 11.42 10.37 -50.37 -12.77 1.06 73 122
## b0 6.87 6.84 1.55 1.48 4.32 9.45 1.00 414 301
## b1 1.88 1.87 0.14 0.14 1.64 2.13 1.01 402 276
## s 2.70 2.81 1.15 1.21 0.77 4.51 1.03 120 38
## sx 1.20 1.21 0.70 0.83 0.19 2.39 1.04 101 115
## x0[1] 4.66 4.66 0.79 0.68 3.32 5.92 1.02 905 1262
## x0[2] 10.62 10.74 0.91 0.79 8.99 12.00 1.01 618 1264
## x0[3] 0.33 0.54 1.24 1.25 -1.97 1.89 1.01 287 381
## x0[4] -0.98 -1.22 1.22 1.27 -2.58 1.14 1.02 197 502
## x0[5] 11.42 11.56 0.95 0.90 9.75 12.76 1.02 361 693
## x0[6] 17.71 17.37 1.51 1.55 15.90 20.72 1.02 165 83
## x0[7] 7.93 7.72 1.16 1.23 6.37 10.07 1.02 200 152
## x0[8] 14.69 14.74 0.84 0.70 13.21 15.99 1.02 841 1998
## x0[9] 1.83 1.99 0.97 0.82 0.05 3.11 1.01 433 453
## x0[10] 14.20 13.94 1.32 1.38 12.54 16.67 1.02 169 61
## x0[11] 12.09 12.07 0.83 0.69 10.75 13.49 1.01 582 740
## x0[12] 2.49 2.30 1.06 1.05 1.05 4.35 1.02 221 416
## x0[13] 10.01 10.15 0.91 0.84 8.42 11.35 1.01 541 1165
## x0[14] 21.44 21.18 1.12 0.90 19.99 23.70 1.01 354 149
## x0[15] 0.69 0.76 0.88 0.65 -0.86 2.05 1.01 990 598
## x0[16] 7.29 7.44 0.92 0.85 5.65 8.59 1.01 415 615
## x0[17] 9.50 9.64 0.89 0.82 7.95 10.80 1.01 511 1262
## x0[18] 5.89 5.94 0.81 0.63 4.46 7.16 1.01 1511 813
## x0[19] 18.43 18.35 0.92 0.71 16.99 20.15 1.01 518 223
## x0[20] 13.56 13.75 1.18 1.30 11.49 15.09 1.02 228 505
## m0[1] 15.63 15.66 1.54 1.39 13.14 18.17 1.00 1618 2603
## m0[2] 26.78 26.81 1.76 1.75 23.99 29.59 1.01 523 2038
## m0[3] 7.54 7.50 2.27 2.66 4.12 11.19 1.01 321 2267
## m0[4] 5.05 5.18 2.50 2.76 0.80 8.66 1.02 183 283
## m0[5] 28.27 28.34 1.84 2.02 25.34 31.14 1.02 295 1981
## m0[6] 40.02 39.85 2.64 3.15 35.96 44.06 1.02 209 1656
## m0[7] 21.73 21.60 2.15 2.48 18.41 25.10 1.02 218 1738
## m0[8] 34.38 34.36 1.66 1.54 31.68 37.13 1.00 1679 2441
## m0[9] 10.33 10.26 1.81 1.87 7.55 13.29 1.01 529 2619
## m0[10] 33.45 33.31 2.34 2.76 29.93 37.09 1.02 232 1933
## m0[11] 29.53 29.56 1.54 1.32 27.00 32.01 1.01 4331 2220
## m0[12] 11.55 11.60 2.08 2.29 8.14 14.76 1.02 214 301
## m0[13] 25.64 25.69 1.75 1.79 22.86 28.45 1.01 484 1835
## m0[14] 47.00 47.15 2.00 1.92 43.59 50.08 1.00 985 2790
## m0[15] 8.20 8.15 1.68 1.47 5.52 11.03 1.01 1428 797
## m0[16] 20.55 20.61 1.74 1.81 17.69 23.32 1.02 355 1589
## m0[17] 24.69 24.69 1.71 1.67 21.93 27.44 1.01 538 2261
## m0[18] 17.92 17.91 1.53 1.37 15.48 20.45 1.00 1885 2411
## m0[19] 41.38 41.41 1.75 1.60 38.48 44.19 1.00 2872 2521
## m0[20] 32.28 32.41 2.34 2.78 28.65 35.92 1.02 209 1863
## y0[1] 15.54 15.57 3.25 2.78 10.22 20.88 1.01 3309 3305
## y0[2] 26.74 26.46 3.43 2.83 21.54 32.71 1.00 2157 3056
## y0[3] 7.58 7.05 3.69 3.42 2.30 14.27 1.01 789 2073
## y0[4] 5.07 5.60 3.85 3.60 -1.74 10.55 1.00 525 1308
## y0[5] 28.32 27.99 3.53 3.16 23.11 34.50 1.00 1031 2844
## y0[6] 39.99 40.48 3.95 3.98 32.78 45.51 1.01 477 3051
## y0[7] 21.71 22.14 3.67 3.43 15.16 26.97 1.01 611 2353
## y0[8] 34.35 34.20 3.32 2.85 29.06 40.02 1.01 3863 2907
## y0[9] 10.37 9.95 3.46 3.01 5.15 16.45 1.01 1859 2480
## y0[10] 33.47 34.01 3.76 3.54 26.82 38.68 1.01 557 2713
## y0[11] 29.55 29.55 3.28 2.75 24.19 35.05 1.01 4263 3552
## y0[12] 11.50 11.92 3.58 3.31 5.20 16.58 1.01 674 2425
## y0[13] 25.61 25.24 3.45 2.98 20.36 31.68 1.00 2123 3397
## y0[14] 46.92 47.25 3.60 3.09 40.74 52.38 1.00 2417 2529
## y0[15] 8.21 8.10 3.45 2.83 2.61 14.05 1.01 3410 2279
## y0[16] 20.52 20.21 3.42 3.05 15.33 26.39 1.00 1461 2592
## y0[17] 24.75 24.41 3.33 2.86 19.81 30.67 1.00 1712 2853
## y0[18] 17.94 17.80 3.33 2.82 12.48 23.55 1.01 3751 3456
## y0[19] 41.35 41.42 3.45 2.79 35.59 47.10 1.00 3917 2795
## y0[20] 32.24 31.67 3.78 3.48 26.98 39.13 1.01 558 1812
## m1[1] 2.98 2.97 1.80 1.73 -0.01 5.95 1.01 409 308
## m1[2] 7.75 7.71 1.50 1.43 5.29 10.24 1.00 414 302
## m1[3] 12.51 12.50 1.23 1.22 10.52 14.58 1.00 417 325
## m1[4] 17.28 17.29 1.02 1.00 15.69 19.01 1.01 421 494
## m1[5] 22.04 22.06 0.91 0.86 20.51 23.52 1.02 412 115
## m1[6] 26.81 26.83 0.94 0.88 25.08 28.35 1.02 447 56
## m1[7] 31.57 31.61 1.09 1.04 29.62 33.34 1.01 418 60
## m1[8] 36.34 36.36 1.33 1.28 33.96 38.54 1.01 390 86
## m1[9] 41.10 41.09 1.62 1.58 38.37 43.82 1.01 386 186
## m1[10] 45.87 45.85 1.93 1.89 42.74 49.12 1.01 385 450
## x1[1] -2.10 -2.09 1.36 0.90 -4.41 0.16 1.02 4210 2527
## x1[2] 0.51 0.50 1.34 0.92 -1.71 2.70 1.02 4279 1693
## x1[3] 2.98 3.00 1.39 0.92 0.69 5.29 1.02 3679 1846
## x1[4] 5.56 5.55 1.38 0.90 3.32 7.87 1.02 3667 2424
## x1[5] 8.09 8.09 1.40 0.92 5.78 10.47 1.02 4233 2071
## x1[6] 10.62 10.63 1.36 0.88 8.34 12.79 1.02 3558 1978
## x1[7] 13.14 13.18 1.39 0.93 10.76 15.38 1.02 3996 2192
## x1[8] 15.72 15.71 1.39 0.94 13.49 18.08 1.02 4210 1766
## x1[9] 18.26 18.23 1.38 0.92 16.07 20.66 1.02 3677 2009
## x1[10] 20.76 20.79 1.38 0.90 18.36 23.01 1.02 3681 1907
## y1[1] 3.08 2.93 3.53 3.15 -2.51 8.91 1.00 1225 1685
## y1[2] 7.75 7.64 3.26 2.82 2.28 13.11 1.00 1422 2882
## y1[3] 12.50 12.40 3.18 2.74 7.29 17.75 1.00 1989 2478
## y1[4] 17.30 17.27 3.15 2.62 12.22 22.47 1.01 3246 1710
## y1[5] 22.02 22.00 3.05 2.62 16.99 27.12 1.00 3624 3149
## y1[6] 26.75 26.72 3.14 2.67 21.68 31.92 1.01 2488 2889
## y1[7] 31.58 31.55 3.17 2.71 26.30 36.74 1.00 2193 2890
## y1[8] 36.37 36.33 3.21 2.91 31.33 41.69 1.00 2054 3064
## y1[9] 41.17 41.09 3.29 3.05 36.04 46.66 1.00 1656 2741
## y1[10] 45.93 45.87 3.61 3.29 40.08 51.78 1.00 999 2168
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -1331.69
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -32.5887 0.00192501 0.000150369 1 1 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -32.59
## b0 7.16
## b1 1.63
## s 3.09
## m0[1] 12.05
## m0[2] 17.06
## m0[3] 10.23
## m0[4] 14.16
## m0[5] 9.01
## m0[6] 18.58
## m0[7] 8.10
## m0[8] 13.12
## m0[9] 17.37
## m0[10] 14.72
## m0[11] 8.89
## m0[12] 16.11
## m0[13] 14.70
## m0[14] 20.16
## m0[15] 15.42
## m0[16] 17.59
## m0[17] 14.98
## m0[18] 17.09
## m0[19] 13.41
## m0[20] 12.08
## y0[1] 17.72
## y0[2] 16.53
## y0[3] 9.67
## y0[4] 14.59
## y0[5] 10.35
## y0[6] 18.42
## y0[7] 7.84
## y0[8] 17.06
## y0[9] 17.46
## y0[10] 18.80
## y0[11] 5.35
## y0[12] 17.40
## y0[13] 12.22
## y0[14] 15.22
## y0[15] 14.65
## y0[16] 16.45
## y0[17] 9.87
## y0[18] 18.11
## y0[19] 13.72
## y0[20] 11.10
## m1[1] 8.10
## m1[2] 9.44
## m1[3] 10.78
## m1[4] 12.12
## m1[5] 13.46
## m1[6] 14.80
## m1[7] 16.14
## m1[8] 17.48
## m1[9] 18.82
## m1[10] 20.16
## y1[1] 5.52
## y1[2] 10.88
## y1[3] 12.94
## y1[4] 18.52
## y1[5] 15.21
## y1[6] 15.75
## y1[7] 10.59
## y1[8] 23.25
## y1[9] 20.31
## y1[10] 22.45
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.06 -32.71 1.33 1.12 -35.66 -31.62 1.02 384 618
## b0 7.22 7.20 1.86 1.73 4.21 10.35 1.02 251 135
## b1 1.61 1.63 0.39 0.37 0.95 2.23 1.01 309 385
## s 3.50 3.42 0.65 0.60 2.61 4.66 1.00 1243 1508
## m0[1] 12.07 12.08 0.96 0.91 10.46 13.66 1.02 325 306
## m0[2] 17.01 17.01 1.04 0.98 15.26 18.73 1.00 1831 1364
## m0[3] 10.27 10.27 1.24 1.15 8.23 12.39 1.02 262 137
## m0[4] 14.15 14.16 0.80 0.75 12.80 15.46 1.00 871 887
## m0[5] 9.06 9.05 1.48 1.36 6.67 11.60 1.02 252 109
## m0[6] 18.52 18.54 1.30 1.26 16.34 20.60 1.00 1004 1270
## m0[7] 8.15 8.14 1.66 1.53 5.46 10.94 1.02 250 121
## m0[8] 13.12 13.14 0.85 0.79 11.68 14.49 1.01 456 443
## m0[9] 17.32 17.33 1.09 1.03 15.50 19.07 1.00 1676 1345
## m0[10] 14.70 14.70 0.81 0.74 13.35 16.04 1.00 1342 1026
## m0[11] 8.93 8.92 1.50 1.38 6.51 11.51 1.02 252 107
## m0[12] 16.08 16.09 0.91 0.84 14.54 17.59 1.00 2234 1403
## m0[13] 14.69 14.69 0.81 0.74 13.34 16.02 1.00 1325 1026
## m0[14] 20.09 20.09 1.61 1.60 17.42 22.62 1.00 697 967
## m0[15] 15.39 15.38 0.85 0.78 13.99 16.79 1.00 2124 1323
## m0[16] 17.54 17.56 1.12 1.07 15.64 19.34 1.00 1504 1343
## m0[17] 14.97 14.96 0.82 0.77 13.60 16.31 1.00 1676 1151
## m0[18] 17.05 17.05 1.04 0.99 15.28 18.77 1.00 1814 1364
## m0[19] 13.41 13.43 0.83 0.78 11.98 14.75 1.01 528 577
## m0[20] 12.10 12.11 0.95 0.90 10.50 13.68 1.02 327 306
## y0[1] 11.96 11.99 3.60 3.48 5.90 17.92 1.00 1107 1363
## y0[2] 16.92 16.93 3.69 3.65 11.05 22.87 1.00 1711 1926
## y0[3] 10.35 10.42 3.69 3.64 4.35 16.44 1.00 1502 1802
## y0[4] 14.17 14.18 3.71 3.53 8.06 20.14 1.00 1886 2011
## y0[5] 9.18 9.16 3.86 3.69 3.06 15.42 1.00 1039 1682
## y0[6] 18.55 18.67 3.80 3.64 12.33 24.85 1.00 1895 2014
## y0[7] 8.12 8.18 4.01 3.85 1.53 14.90 1.00 996 1239
## y0[8] 12.98 13.01 3.68 3.56 7.08 18.97 1.00 1768 1915
## y0[9] 17.35 17.42 3.86 3.62 11.24 23.67 1.00 1656 1955
## y0[10] 14.50 14.48 3.63 3.47 8.58 20.57 1.00 2033 1961
## y0[11] 8.97 8.98 3.75 3.64 2.70 15.03 1.00 981 1730
## y0[12] 16.25 16.21 3.71 3.42 10.19 22.38 1.00 2224 2008
## y0[13] 14.79 14.77 3.71 3.54 8.63 20.84 1.00 1917 2055
## y0[14] 20.16 20.08 3.99 3.82 13.76 26.88 1.00 1485 1699
## y0[15] 15.35 15.33 3.64 3.52 9.28 21.06 1.00 1972 1839
## y0[16] 17.52 17.50 3.91 3.75 11.20 23.88 1.00 1943 1893
## y0[17] 14.79 14.75 3.67 3.51 8.70 20.82 1.00 1827 1890
## y0[18] 17.08 17.03 3.83 3.71 10.88 23.25 1.00 2067 1971
## y0[19] 13.40 13.35 3.70 3.63 7.31 19.80 1.00 1974 1904
## y0[20] 12.00 12.00 3.72 3.51 5.96 17.98 1.00 1526 1745
## m1[1] 8.15 8.14 1.66 1.53 5.46 10.94 1.02 250 121
## m1[2] 9.48 9.48 1.39 1.28 7.23 11.86 1.02 254 108
## m1[3] 10.81 10.83 1.15 1.08 8.91 12.76 1.02 272 171
## m1[4] 12.13 12.14 0.95 0.90 10.54 13.70 1.02 329 306
## m1[5] 13.46 13.48 0.82 0.78 12.03 14.80 1.01 542 571
## m1[6] 14.78 14.78 0.81 0.75 13.42 16.12 1.00 1437 1104
## m1[7] 16.11 16.11 0.92 0.84 14.56 17.62 1.00 2226 1403
## m1[8] 17.43 17.44 1.11 1.05 15.57 19.22 1.00 1595 1344
## m1[9] 18.76 18.78 1.34 1.31 16.52 20.91 1.00 933 1132
## m1[10] 20.09 20.09 1.61 1.60 17.42 22.62 1.00 697 967
## y1[1] 8.08 8.09 3.98 3.77 1.64 14.43 1.00 762 1178
## y1[2] 9.53 9.57 3.79 3.74 3.32 15.66 1.00 1078 1292
## y1[3] 10.74 10.75 3.70 3.68 4.73 16.72 1.00 1506 1788
## y1[4] 12.12 12.11 3.70 3.53 5.99 18.27 1.00 1851 1877
## y1[5] 13.43 13.47 3.63 3.45 7.32 19.27 1.00 2014 1694
## y1[6] 14.81 14.81 3.51 3.42 8.87 20.50 1.00 1988 1673
## y1[7] 16.15 16.20 3.63 3.53 10.15 22.01 1.00 1775 1835
## y1[8] 17.61 17.59 3.86 3.74 11.52 23.83 1.00 1967 1866
## y1[9] 18.74 18.75 3.68 3.45 12.76 24.59 1.00 1828 1674
## y1[10] 20.05 20.01 3.89 3.80 13.60 26.38 1.00 1448 1646
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -110.328
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -10.1859 6.68021e-05 0.000461082 1 1 32
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -10.19
## b0 5.27
## b1 1.89
## s 0.54
## m0[1] 10.93
## m0[2] 16.71
## m0[3] 8.83
## m0[4] 13.36
## m0[5] 7.41
## m0[6] 18.47
## m0[7] 6.36
## m0[8] 12.16
## m0[9] 17.07
## m0[10] 14.01
## m0[11] 7.27
## m0[12] 15.62
## m0[13] 13.99
## m0[14] 20.30
## m0[15] 14.82
## m0[16] 17.33
## m0[17] 14.32
## m0[18] 16.75
## m0[19] 12.50
## m0[20] 10.97
## y0[1] 11.09
## y0[2] 17.45
## y0[3] 8.81
## y0[4] 13.56
## y0[5] 7.19
## y0[6] 17.92
## y0[7] 7.66
## y0[8] 10.82
## y0[9] 19.23
## y0[10] 13.67
## y0[11] 7.02
## y0[12] 16.29
## y0[13] 18.61
## y0[14] 20.06
## y0[15] 14.97
## y0[16] 17.67
## y0[17] 13.19
## y0[18] 17.10
## y0[19] 12.77
## y0[20] 10.52
## m1[1] 6.36
## m1[2] 7.91
## m1[3] 9.46
## m1[4] 11.01
## m1[5] 12.55
## m1[6] 14.10
## m1[7] 15.65
## m1[8] 17.20
## m1[9] 18.75
## m1[10] 20.30
## y1[1] 6.35
## y1[2] 8.88
## y1[3] 10.17
## y1[4] 10.00
## y1[5] 12.39
## y1[6] 15.03
## y1[7] 15.20
## y1[8] 16.75
## y1[9] 18.97
## y1[10] 20.18
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.37 -12.04 1.25 0.99 -14.86 -10.98 1.01 602 1090
## b0 5.30 5.26 0.50 0.51 4.55 6.20 1.01 474 589
## b1 1.89 1.89 0.11 0.11 1.71 2.06 1.01 588 702
## s 0.64 0.62 0.19 0.18 0.38 1.01 1.00 760 837
## m0[1] 10.97 10.94 0.26 0.24 10.58 11.41 1.00 636 777
## m0[2] 16.76 16.74 0.30 0.29 16.29 17.26 1.00 2113 1268
## m0[3] 8.86 8.83 0.34 0.33 8.36 9.45 1.01 501 613
## m0[4] 13.40 13.39 0.23 0.22 13.06 13.78 1.00 1520 1431
## m0[5] 7.45 7.41 0.40 0.40 6.85 8.14 1.01 479 568
## m0[6] 18.52 18.49 0.38 0.36 17.95 19.16 1.00 1510 1094
## m0[7] 6.39 6.35 0.45 0.45 5.72 7.19 1.01 475 571
## m0[8] 12.20 12.18 0.23 0.22 11.84 12.61 1.00 903 1011
## m0[9] 17.12 17.09 0.32 0.30 16.64 17.66 1.00 2023 1251
## m0[10] 14.05 14.04 0.23 0.22 13.69 14.44 1.00 2139 1435
## m0[11] 7.30 7.26 0.41 0.41 6.70 8.01 1.01 478 568
## m0[12] 15.67 15.65 0.27 0.25 15.26 16.11 1.00 2315 1604
## m0[13] 14.03 14.02 0.23 0.22 13.68 14.42 1.00 2131 1456
## m0[14] 20.35 20.31 0.47 0.45 19.64 21.16 1.00 1181 928
## m0[15] 14.86 14.85 0.24 0.23 14.48 15.28 1.00 2326 1578
## m0[16] 17.38 17.35 0.33 0.31 16.89 17.93 1.00 1950 1204
## m0[17] 14.36 14.35 0.24 0.23 14.00 14.76 1.00 2232 1458
## m0[18] 16.80 16.78 0.30 0.29 16.33 17.31 1.00 2105 1268
## m0[19] 12.54 12.52 0.23 0.22 12.19 12.93 1.00 1044 1169
## m0[20] 11.00 10.98 0.26 0.24 10.61 11.44 1.00 640 777
## y0[1] 10.53 10.96 27.11 0.93 7.63 15.48 1.00 1842 1866
## y0[2] 17.57 16.78 33.14 0.96 12.93 21.03 1.00 1637 1973
## y0[3] -2.84 8.86 521.28 1.03 5.45 13.28 1.00 1852 1930
## y0[4] 14.00 13.43 34.84 0.91 9.73 17.96 1.00 2046 1890
## y0[5] 7.07 7.41 18.61 1.13 3.34 11.94 1.00 1655 1723
## y0[6] 17.99 18.53 43.71 1.02 15.18 22.07 1.00 1923 1915
## y0[7] 6.69 6.35 38.70 1.12 2.14 11.44 1.00 1910 1725
## y0[8] 10.29 12.14 96.11 0.93 8.12 16.44 1.00 1808 1822
## y0[9] 20.62 17.08 203.42 0.97 13.22 20.90 1.00 1958 1932
## y0[10] 14.00 14.04 24.09 0.93 10.41 17.75 1.00 1980 1972
## y0[11] 4.83 7.27 199.45 1.06 3.03 11.69 1.00 1590 1760
## y0[12] 15.34 15.67 24.00 1.02 11.86 20.25 1.00 2164 1890
## y0[13] 15.15 14.02 43.40 0.98 10.06 18.25 1.00 1994 1926
## y0[14] 20.33 20.35 21.39 1.13 16.15 24.46 1.00 1808 1971
## y0[15] 14.02 14.84 13.89 0.93 10.22 18.37 1.00 2045 1897
## y0[16] 17.90 17.36 19.60 1.06 13.24 21.22 1.00 1874 1898
## y0[17] 14.59 14.32 83.19 0.97 10.07 19.16 1.00 2031 1974
## y0[18] 16.82 16.81 16.24 1.01 12.77 21.34 1.00 2089 1883
## y0[19] 12.96 12.52 33.13 0.97 8.49 16.42 1.00 2019 1932
## y0[20] 6.24 11.03 203.23 0.99 7.42 15.24 1.00 1930 1931
## m1[1] 6.39 6.35 0.45 0.45 5.72 7.19 1.01 475 571
## m1[2] 7.94 7.91 0.38 0.37 7.38 8.59 1.01 483 582
## m1[3] 9.49 9.46 0.31 0.30 9.03 10.04 1.01 523 691
## m1[4] 11.04 11.02 0.26 0.24 10.66 11.48 1.00 645 777
## m1[5] 12.60 12.58 0.23 0.22 12.25 12.99 1.00 1074 1197
## m1[6] 14.15 14.14 0.23 0.22 13.79 14.54 1.00 2167 1458
## m1[7] 15.70 15.68 0.27 0.25 15.29 16.14 1.00 2311 1604
## m1[8] 17.25 17.23 0.32 0.31 16.77 17.80 1.00 1983 1251
## m1[9] 18.80 18.78 0.39 0.37 18.21 19.47 1.00 1440 1107
## m1[10] 20.35 20.31 0.47 0.45 19.64 21.16 1.00 1181 928
## y1[1] 6.61 6.39 7.80 1.11 2.41 11.02 1.00 1514 1706
## y1[2] 7.66 7.90 9.53 0.97 4.15 11.28 1.00 1827 1926
## y1[3] 9.61 9.51 33.61 0.96 5.58 13.27 1.00 1936 2151
## y1[4] 10.87 11.05 22.07 0.97 7.59 14.77 1.00 1735 1888
## y1[5] 12.51 12.57 8.93 0.95 8.44 16.11 1.00 1928 1984
## y1[6] 14.06 14.13 16.35 0.96 10.91 18.08 1.00 2055 1899
## y1[7] 15.67 15.68 17.00 0.96 11.35 20.10 1.00 1891 2059
## y1[8] 17.19 17.23 12.93 1.10 12.43 21.75 1.00 1908 1984
## y1[9] 19.32 18.82 27.96 1.00 14.71 22.52 1.00 1857 1866
## y1[10] 21.78 20.35 56.80 1.05 16.13 24.84 1.00 1489 1955
censored data
y i=1-N, y~N(m,s)
acutualy ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11.stan
objective variable is censored
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitiondata=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -6470.48
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 -13.5528 0.00361578 0.000438694 1 1 40
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.55
## b0 17.65
## b1 1.11
## s 2.73
## m0[1] 21.39
## m0[2] 21.89
## m0[3] 25.42
## m0[4] 24.32
## m0[5] 23.65
## m0[6] 23.82
## m0[7] 22.15
## m0[8] 21.65
## m0[9] 20.24
## y0[1] 18.06
## y0[2] 21.57
## y0[3] 26.36
## y0[4] 28.30
## y0[5] 23.82
## y0[6] 28.16
## y0[7] 21.27
## y0[8] 23.55
## y0[9] 17.85
## m1[1] 17.92
## m1[2] 18.99
## m1[3] 20.07
## m1[4] 21.14
## m1[5] 22.22
## m1[6] 23.29
## m1[7] 24.37
## m1[8] 25.44
## m1[9] 26.52
## m1[10] 27.59
## y1[1] 18.16
## y1[2] 16.36
## y1[3] 20.64
## y1[4] 23.39
## y1[5] 22.29
## y1[6] 20.96
## y1[7] 17.38
## y1[8] 23.67
## y1[9] 24.31
## y1[10] 28.30mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.41 -14.07 1.48 1.22 -17.74 -12.76 1.01 308 295
## b0 17.83 17.95 4.59 4.39 10.67 25.19 1.03 227 156
## b1 1.07 1.04 0.97 0.89 -0.48 2.53 1.03 244 200
## s 3.89 3.57 1.33 1.06 2.36 6.48 1.00 649 615
## m0[1] 21.42 21.45 1.74 1.64 18.52 24.16 1.01 290 194
## m0[2] 21.91 21.96 1.51 1.38 19.42 24.27 1.01 381 402
## m0[3] 25.31 25.28 2.81 2.44 20.83 29.73 1.02 388 415
## m0[4] 24.25 24.25 2.02 1.82 21.06 27.47 1.01 562 478
## m0[5] 23.60 23.61 1.63 1.45 21.05 26.18 1.01 888 710
## m0[6] 23.77 23.76 1.72 1.53 21.06 26.48 1.01 773 613
## m0[7] 22.16 22.21 1.42 1.28 19.80 24.44 1.00 477 557
## m0[8] 21.68 21.72 1.61 1.50 19.03 24.20 1.01 326 299
## m0[9] 20.32 20.39 2.50 2.42 16.20 24.30 1.02 237 146
## y0[1] 21.39 21.50 4.63 4.05 13.67 28.75 1.00 1098 901
## y0[2] 21.65 21.74 4.46 3.95 14.36 28.80 1.00 1567 1703
## y0[3] 25.23 25.25 4.96 4.27 17.23 32.93 1.01 761 523
## y0[4] 24.29 24.30 4.57 3.97 16.80 31.31 1.01 1549 1440
## y0[5] 23.54 23.51 4.60 3.92 16.37 30.68 1.00 2014 1653
## y0[6] 23.87 23.83 4.50 3.87 16.67 30.99 1.00 1863 1892
## y0[7] 22.12 22.17 4.32 3.79 15.25 28.81 1.00 1362 1141
## y0[8] 21.58 21.61 4.43 3.89 14.33 28.53 1.00 1646 1476
## y0[9] 20.25 20.13 4.76 4.37 12.91 27.89 1.00 806 723
## m1[1] 18.09 18.19 4.36 4.21 11.30 25.11 1.03 226 156
## m1[2] 19.12 19.16 3.48 3.34 13.62 24.73 1.03 228 157
## m1[3] 20.15 20.20 2.63 2.55 15.84 24.32 1.02 235 155
## m1[4] 21.19 21.25 1.89 1.80 18.06 24.10 1.02 267 126
## m1[5] 22.22 22.25 1.41 1.26 19.88 24.47 1.00 511 607
## m1[6] 23.26 23.25 1.48 1.30 20.93 25.59 1.00 1146 872
## m1[7] 24.29 24.28 2.05 1.84 21.06 27.57 1.01 550 469
## m1[8] 25.33 25.30 2.82 2.46 20.82 29.80 1.02 386 415
## m1[9] 26.36 26.32 3.68 3.25 20.47 32.06 1.02 330 347
## m1[10] 27.40 27.34 4.57 4.07 20.05 34.52 1.02 304 278
## y1[1] 18.04 18.05 5.92 5.35 8.43 27.41 1.02 382 388
## y1[2] 19.10 19.15 5.46 4.97 10.15 27.98 1.01 459 433
## y1[3] 20.01 20.17 4.98 4.28 11.49 27.83 1.01 645 915
## y1[4] 21.11 21.04 4.63 4.05 13.69 28.35 1.01 1032 1378
## y1[5] 22.37 22.46 4.29 3.77 15.21 29.10 1.00 1955 1771
## y1[6] 23.21 23.24 4.33 3.87 16.46 29.97 1.00 2048 1796
## y1[7] 24.15 24.20 4.52 3.78 16.81 31.51 1.00 1846 1274
## y1[8] 25.38 25.49 4.98 4.47 17.19 33.30 1.01 1012 875
## y1[9] 26.42 26.59 5.64 4.89 17.31 35.16 1.01 660 705
## y1[10] 27.54 27.41 6.22 5.35 17.77 37.71 1.01 494 787
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -68.8876
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 25 -18.2624 0.00137193 0.000255376 1 1 30
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -18.26
## b0 9.66
## b1 2.90
## s 3.32
## m0[1] 19.47
## m0[2] 20.79
## m0[3] 30.04
## m0[4] 27.17
## m0[5] 25.39
## m0[6] 25.85
## m0[7] 21.48
## m0[8] 20.15
## m0[9] 16.46
## y0[1] 17.57
## y0[2] 17.99
## y0[3] 28.02
## y0[4] 29.44
## y0[5] 23.86
## y0[6] 29.86
## y0[7] 22.77
## y0[8] 22.97
## y0[9] 18.34
## m1[1] 10.37
## m1[2] 13.19
## m1[3] 16.01
## m1[4] 18.83
## m1[5] 21.64
## m1[6] 24.46
## m1[7] 27.28
## m1[8] 30.10
## m1[9] 32.92
## m1[10] 35.74
## y1[1] 10.95
## y1[2] 10.25
## y1[3] 16.16
## y1[4] 22.28
## y1[5] 22.72
## y1[6] 24.66
## y1[7] 32.29
## y1[8] 29.59
## y1[9] 29.18
## y1[10] 33.12mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -18.82 -18.39 1.52 1.20 -21.89 -17.25 1.00 355 466
## b0 6.89 7.75 4.47 3.50 -1.36 12.12 1.01 303 306
## b1 3.52 3.33 0.91 0.69 2.47 5.14 1.01 317 358
## s 4.77 4.38 1.84 1.28 2.84 7.99 1.00 512 567
## m0[1] 18.79 19.01 1.87 1.50 15.45 21.37 1.01 417 472
## m0[2] 20.40 20.52 1.64 1.31 17.57 22.76 1.01 519 592
## m0[3] 31.63 31.19 2.60 2.03 28.52 36.15 1.00 589 664
## m0[4] 28.15 27.87 1.93 1.57 25.75 31.44 1.00 961 778
## m0[5] 25.99 25.81 1.62 1.36 23.83 28.71 1.00 1433 945
## m0[6] 26.54 26.33 1.69 1.42 24.34 29.37 1.00 1321 938
## m0[7] 21.23 21.31 1.56 1.27 18.59 23.48 1.01 610 647
## m0[8] 19.62 19.78 1.74 1.42 16.58 22.07 1.01 459 531
## m0[9] 15.14 15.56 2.56 2.08 10.49 18.36 1.01 334 360
## y0[1] 18.82 18.92 5.44 4.60 9.47 27.51 1.00 1672 1002
## y0[2] 20.36 20.37 5.31 4.66 11.87 28.97 1.00 1914 1464
## y0[3] 31.80 31.34 5.93 5.15 23.39 41.99 1.00 1508 1213
## y0[4] 28.02 27.68 5.45 4.52 19.94 36.66 1.00 1881 1875
## y0[5] 25.87 25.83 5.26 4.49 17.84 34.12 1.00 2092 1610
## y0[6] 26.34 26.17 5.28 4.62 18.19 34.91 1.00 1763 1500
## y0[7] 21.37 21.30 5.44 4.56 13.03 29.60 1.00 1846 1633
## y0[8] 19.79 20.04 5.19 4.55 11.64 27.80 1.00 1628 1157
## y0[9] 15.13 15.60 5.53 4.94 4.95 23.11 1.01 980 983
## m1[1] 7.75 8.58 4.26 3.32 0.04 12.78 1.01 304 307
## m1[2] 11.17 11.82 3.45 2.68 4.92 15.31 1.01 310 328
## m1[3] 14.59 15.01 2.67 2.17 9.68 17.92 1.01 328 364
## m1[4] 18.01 18.25 2.00 1.64 14.39 20.73 1.01 388 417
## m1[5] 21.44 21.51 1.54 1.25 18.84 23.68 1.01 639 657
## m1[6] 24.86 24.73 1.52 1.29 22.71 27.36 1.00 1439 867
## m1[7] 28.28 27.99 1.95 1.58 25.87 31.61 1.00 937 730
## m1[8] 31.70 31.25 2.62 2.04 28.57 36.23 1.00 585 664
## m1[9] 35.12 34.51 3.38 2.63 31.10 41.11 1.01 466 488
## m1[10] 38.54 37.72 4.19 3.27 33.58 45.84 1.01 415 472
## y1[1] 7.74 8.32 6.43 5.59 -3.67 16.94 1.00 799 650
## y1[2] 11.05 11.65 6.17 5.23 -0.02 19.76 1.00 690 825
## y1[3] 14.50 14.89 5.72 4.97 4.60 23.24 1.00 968 1022
## y1[4] 17.96 18.26 5.66 4.69 8.44 26.32 1.00 1226 1235
## y1[5] 21.56 21.59 5.32 4.51 13.26 29.97 1.00 1332 1572
## y1[6] 24.98 24.83 5.35 4.37 16.82 34.16 1.00 1656 1233
## y1[7] 28.15 27.85 5.43 4.85 20.00 37.13 1.00 1814 1326
## y1[8] 31.70 31.19 5.74 4.69 23.53 42.43 1.00 1351 1023
## y1[9] 35.00 34.47 6.31 5.31 26.37 45.51 1.00 901 932
## y1[10] 38.47 37.82 6.32 5.38 29.82 49.32 1.00 885 915
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
appendix: bimodal distribution, mixed normal distribution
stan
data {
int<lower=0> N;
real y[N];
}
parameters {
real<lower=0, upper=1> theta; //ratio of mix
real mu1;
real mu2;
real<lower=0> sigma1;
real<lower=0> sigma2;
}
model {
for (n in 1:N) {
target += log_mix(theta,
normal_lpdf(y[n] | mu1, sigma1),
normal_lpdf(y[n] | mu2, sigma2));
}
}EM algorithm
y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -242 100 6 -511 -512
##
## Clustering table:
## 1 2 3
## 35 34 31rst$parameters## $pro
## [1] 0.356 0.334 0.310
##
## $mean
## 1 2 3
## -4.8267 0.0605 5.5018
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 0.841plot(rst)
